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Abstract Numerical calculation and analysis of extremely high-lying energy spectra, containing 
thousands of levels with sequential quantum number up to 62,000 per symmetry class, of a 
generic chaotic 3D quantum billiard is reported. The shape of the billiard is given by a simple 
and smooth de formation of a unit sphere which gives rise to (almost) fully chaotic classical 
dynamics. We present an analysis of (i) quantum length spectrum whose smooth part agrees 
with the 3D Weyl formula and whose oscillatory part is peaked around the periods of classical 
periodic orbits, (ii) nearest neighbor level spacing distribution and (iii) number variance. 

Although the chaotic classical dynamics quickly and uniformly explores almost entire energy 
shell, while the measure of the regular part of phase space is insignificantly small, we find small 
but significant deviations from GOE statistics which are explained in terms of localization of 
eigenf unctions onto lower dimensional classically invariant manifolds. 
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In recent years a vast amount of numerical work has been invested in understanding the con- 
nection between classical chaotic dynamics and statistical properties of quantum energy spectra 
of autonomous Hamiltonian systems (see e.g. j ij and references therein). Theoretical results 
relate the statistical properties of quantum energy levels either with the classical periodic orbits 
or with the eigenvalues of the statistical ensembles of random matrices [|| . In the later sense the 
strongest result is numertically and theoretically supported conjecture by Bohigas, Giannoni and 
Schmit |J which identifies the asymptotic statistical properties of quantal energy spectra of clas- 
sically fully chaotic systems with those of ensembles of Gaussian orthogonal/unitary/symplectic 
(GOE/GUE/GSE) random matrices; fine scale deviations from expected asymptotic, say GOE, 
statistics are related to localization — non-uniformity of quantum eigenstates on the energy 
surface. So far, people have mostly studied systems with two degrees of freedom which are 
numerically much more accessible than higher dimensional systems, but which are also rather 
special since in 2D systems, unlike in higher dimensional systems, invariant tori (or cantori - 
partial barriers) can geometrically divide the energy surfaces. However, numerics in quantum 
3D systems is much more computer time consuming than in 2D systems, so it is very difficult 
to calculate highly excited eigenstates or long spectral stretches of chaotic 3D systems. Due to 
their energy-scaling property the cleanest systems to work with are certainly the 3D billiards — 
free particles moving freely in a bounded 3D domain and bouncing off its boundary elastically. 
Recently, Primack and Smilansky 0] have calculated the first two thousand eigenstates of a 
3D Sinai billiard, the free particle in the region inside a cube and outside a smaller concentric 
sphere. Although the 3D Sinai billiard is classically fully chaotic — ergodic, it is nongeneric in 
a sense that it possesses 2D and 3D invariant manifolds of non-isolated "bouncing ball" periodic 
orbits which greatly affect the properties of quantal spectra and produce significant deviations 
from GOE statistics for finite spectral samples as explained in 

The first aim of this paper is to consider spectral statistics of a generic classically chaotic 
3D system. We have chosen a 3D billiard with a smooth C°° boundary. Since no such system 
is known to be rigorously ergodic we have defined a two-parameter family of 3D billiards whose 
shapes are given by simplest smooth deformations of a sphere: The radial distance rs(n) from 
the origin to the boundary is defined to be the following function of the direction n,n 2 = 1, 

re(^) = 1 + a(n x + riy + n z ) + bri^riyn 2 (1) 

which contains the two lowest order terms which preserve the cubic symmetry (the first two 
fully symmetric type (a) cubic harmonics after Von Lage and Bethe ||). Using efficiently coded 
billiard's classical dynamics we have numerically explored the parameter space (a, b) and found 
generic transition from integrable sphere, for a = b = 0, via KAM scenario to more and more 
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chaotic cases for increasing absolute values of parameters a and b. Finally we have chosen the 
shape (see figure 1): a = — l/5,b = —12/5 for which the billiard has the following proper- 
ties: (i) Almost all orbits of the system are uniformly and strongly chaotic with large average 
maximal liapunov exponent, {\ max ) = 0.54, and sharply peaked distribution of maximal finite- 
time liapunov exponets even for a short time of simulation per orbit (see figure 2). There is 
no diffusion-like transport or partial barriers in phase space, (ii) The relative measure of the 
regular part of phase space pi ~ 1CP 3 is very small and should be almost insignificant for quan- 
tal calculations (see figure 2). (iii) The billiard is convex, (iv) There are periodic orbits (for 
unspecified velocity they are (1+1)D invariant submanifolds) which support all types of sta- 
bility: besides the hyperbolic (including loxodromic) and a small fraction of elliptic orbits the 
system possesses also neutrally stable isolated parabolic orbits which touch the boundary only 
at the points of zero curvature radii, (±ri, 0, 0), (0, ±rj., 0), (0, 0, ±ri), r\ =r^(l,0,0). (v) There 
are (2+2)D invariant submanifolds in phase space (r,p). Besides the trivial symmetry planes 
{x = 0,p x = 0} and {x = y,p x = Py} (and manifolds corresponding to cyclic permutations 
of axes) and the boundary surface {^(r/r) = r,p- V(re(r/r) — r) = 0}, there is a nontrivial 
invariant submanifold {x + y = ±z,p x +p y = ±p 2 } (and its cyclic permutations) which is not a 
symmetry plane. The existence of such nontrivial flat invariant submanifold requires that on the 
curve, which is an intersection of the plane x + y = ±z and the boundary surface, any normal 
to the boundary surface should lie in the plane x + y = ±z. For (|l]) this is true if b = 12a. The 
existence and construction of non-flat (2+2) invariant submanifolds is to the best of author's 
knowledge an open problem. 

The second aim of this paper is to demonstrate the power of the scaling method for quantization 
of billiards, proposed recently by Vergini and Saraceno f|], in three dimensions. We have slightly 
modified their method by using expansions in terms of spherical waves (ftkim( r n) = ji(kr)Yi m (n) 
instead of plane waves exp(ik ■ rn), the former being much closer to billiard's geometry and 
naturally incorporating the evanescent modes which are necessary to make the method arbi- 
trarily accurate. The number of spherical waves necessary to accurately quantize a 3D billiard 
for wavenumbers around k, Nsw(k) ~ (kr max ) 2 , where r max = max{rs(n);n 2 = 1}, is almost 
optimal for shapes close to a sphere. The efficiency of the basis of spherical scaling functions 
{4>kim\ with respect to a 3D billiard with the shape (||) is defined as 

M{k) _ V 
V M max {k) ~ 4vrr3 iaa ,/3 

where M{k) = #{k n < k} ~ Vk 3 /(6ir 2 ) (see Weyl formula below (||)) is the number of levels 
with wavenumber k n less than k and V = (1/3) / d 2 nr^ B {n) is the volume of the billiard, while 
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N max (k) = 2(kr max ) 3 / (9tt) is the number of levels less than k in the smallest enscribed spherical 
billiard which has radius r max . For the spherical billiard efficiency rj is equal to one, whereas in 
order to have efficient method of quantization of a strongly chaotic system we have to compro- 
mise between largest possible efficiency n, on one hand, and largest possible average liapunov 
exponent (X max ) with small measure of regular orbits pi, on the other hand. We believe that our 
shape a = —1/5,6 = —12/5, which has efficiency rj = 0.866, represents in this respect somehow 
an optimal choice. More details on the scaling approach to quantization of 3D billiards can be 
found in a subsequent paper |7]]. 

The billiard is preserved under 48-fold cubic symmetry group, so we have desymmetrized it 
and chosen only the class of singlet states belonging to fully symmetric irrep of a cubic group, 
i.e. we have considered a desymetrized part of 1/48 of a billiard bounded by the symmetry 
planes where von Neuman boundary conditions are imposed and the boundary surface where 
Dirichled boundary conditions are imposed. Desymmetrization also reduces the number of scal- 
ing functions, roughly by a factor of 48, which now become the cubic waves (CW): products of 
spherical Bessel functions and cubic harmonics of type (a) of Von Lage and Bethe [|| (i.e. linear 
combinations of spherical harmonics Y[ m with fixed I which are invariant under the actions of 
the cubic group Oh), Ncw{k) ~ (kr max ) 2 /48. (See |?J for more details on the computation of 
CW.) We have calculated two spectral stretches containing around 7000 consecutive levels each: 
sample I contained 6973 consecutive levels on the interval 99.37 < k < 200.38, and sample II 
contained 7133 consecutive levels on the interval 384.23 < k < 400.46. In a single computer run, 
which required to solve a generalized eigenvalue problem |], [?J with matrices of size Ncw(k), 
and has taken about 6 hours of Convex C3860 CPU time we have obtained around 250 accurate 
consecutive levels. (Few tens comparable computer jobs would be required to determine a single 

level with Heller's plane wave decomposition or boundary integral method.) Maximal error was 
allowed to be 10 -3 of the mean level spacing (MLS) but inaccuracy for most of eigenvalues was 
between 10 -5 MLS and 10 -4 MLS. It seemed that the number of accurate levels N conver g e d(k) 
per diagonalization is indeed proportional with the dimension of matrices Ncw(k) as claimed 
by Vergini and Saraceno Q, 

N converg ed{k) « 0.1Ncw(k) « 0.0016& 2 . 

The scaling method has no problem of missing levels. We have compared the smooth part of 
the level counting function M(E) with the generalization of the popular Weyl formula Q from 
2D to 3D billiards 

M smooth {k) = V k s + A N -A D + 

OTT z 167T 
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Uhr> ^ 2q + , 2a 

where V is the volume, An (Ad) are the total areas of the boundary surfaces where Von Neuman 
(Dirichlet) boundary conditions are imposed. It is assumed that pairs of neighbouring smooth 
parts of the boundary surface with the same (opposite) type of boundary condition, DD or NN 
(DN or ND), intersect under the constant angle (a^) along the curve of length (Z^). C 
is the surface integral of the trace of curvature matrix 



c = I d2s (i + i ] > 

J JXl 1X2 



boundary 

where R\ and R2 are the main curvature radii. In figure 3 we plot the difference N(E) — 
M smoo th{E) for the sample I which locally average to zero so well that one could detect an error 
of 5% in the curvature term C\ Due to formula @ the sample I contains the energy levels from 
1,047th to 8,018th excited state, while the sample II contains energy levels from 54,706th to 
61,842th excited state with uncertainty ±2. 

In order to analyze the relation of the spectra to classical periodic orbits we have calculated 
the quantum length spectrum (slightly modified w.r.t refs. P, BJ]) 



k-2 

D(x) = J dk w(k; k\, £2) cos(kx)d osc (k) (3) 



I 

ki 

where d osc (k) = (d/ dk)(M(k) — M smoo th(k)) is an oscillatory part of the density of states and 

(k 2 -k)(k-k 1 ) 



w(k; k\, h 



2) 



6(k 2 - hf 



is a Welsch window function which reduce oscillations due to final spectral length. The length 
spectrum can be expressed semiclassically |J in terms of a sum of fat delta functions of width 
~ l/(k2 — k\) peaked around the periods of the classical periodic orbits (of a desimetrized 
billiard) with weights which are inversely proportional to the stability exponents. We have cal- 
culated 3692 periodic orbits of a desymmetrized billiard for up to 12 bounces and managed to 
identify most of the peaks in the length spectrum of sample I with the least unstable periodic 
orbits (figure 4a). Few unmatched peaks correspond to either missed or longer orbits. For the 
sample II the resolution 5x = 2n/(k2 — k\) is much worse and the broader peaks can no longer be 
associated with individual periodic orbits but with the families of similar periodic orbits which 
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have similar length and therefore can interfere constructively (figure 4b). 

To check the predictions of random matrix theory, the spectra {k n } were unfolded k n — > e n 
to unit MLS, e n = ftf sm ooth{k n ) ■ First we have studied the nearest neighbor level spacing dis- 
tribution (NNLSD) P(S) of spacings S n — c-n+i — &n and compared it with the predictions of 
GOE (accurate approximation of the exact infinitely dimensional formula for Pgoe(S) fl(| have 
been used rather than the popular 2-dim Wigner surmise). We have analyzed cumulative level 
spacing distribution W(S) = Jq dsP(s) = #{5* n < S}/#{S n } in order to avoid losing infor- 
mation by bining and we have found small but statistically significant deviations of NNLSD 
from Pgoe(S) for both samples I and II (see figure 5). Deviations for the higher sample II, 
which are represented with uniform statistical error, U(Wn(S)) — U(Wgoe(S)) (see figure 5b), 
where U(W) = (2 fir) arccos Vl — W (see [|ll]]), are functionally similar and only slightly smaller 
than deviations for the lower sample I, U(Wi(S)) — U(Wgoe(S)). (So deviations are clearly 
of systematic and not of statistical origin.) We have tried to fit the NNLSD also with the 
phenomenological Brody distribution which successfully describes the so-called fractional 
power law level repulsion ]ll|], Pr(S — > 0) oc S^, typically associated with a localization of 
eigenstates. We have also tried the two component Berry-Robnik NNLSD P pi (S) |0| which 
assumes that the spectrum can be decomposed as statistically independent superposition of 
an uncorrelated level sequence with measure p\ which may correspond to localized or regular 
states and GOE level sequence with measure p2 = 1 — pi which may correspond to extended — 
delocalized chaotic states. As a result of the least square fit (see also figure 5b) we have obtained 
for the sample I: (3i = 0.872, pn = 0.024 and (slightly more towards GOE) for the sample II: 
Pu = 0.892, p\n = 0.019. (The Brody fit is with both samples better than the Berry-Robnik 
fit, but neither of them is statistically significant.) 

The NNLSD depends mostly on the short range correlations between quantum energy lev- 
els. To compare the long range correlations with the predictions of random matrix theory 

we have studied the number variance S 2 (L) = (^(J\f(e + L) — J\f{e) — L) 2 ^, the variance of the 

number of unfolded levels in the interval of length L, where N{e) = #{e n < e} counts un- 
folded levels below e. For the GOE spectra the number variance should increase logarithmically 
Yiqq E (L) ~ -Tj log(2-7rL) whereas for the uncorrelated (Poissonian) sequence of levels the num- 
ber variance is equal to the number itself S 2 ncorr (L) = L. We have again found significant 
deviations from GOE statistics for both spectral samples but deviations for the lower sample 
I were significantly larger than for the higher sample II which is consistent with the validity 
of GOE statistics in the ultimate semiclassical limit. We conclude that localized states consti- 



6 



tute approximately independent uncorrected spectral subsequence, with respect to long-range 
correlations, since the number variance is excellently captured by the ansatz of Seligman and 
Verbaarschot |l4j (see figure 6) 

Y?{L)=p l L + Y? GOE (~p 2 L), (4) 

with the following measures of uncorrelated (localized) subsequence pu = 0.072 and pm = 
0.028, for the samples I and II, respectively. In other words: in the samples I and II there should 
be « 500 and as 200 strongly localized eigenstates (scars), respectively [0]. However, note that 
one should not expect consistency with Berry-Robnik parameter p\ (as determined from non- 
significant fit by the Berry-Robnik distribution ]l3[), since the two spectral subsequences are 
expected to be short-range correlated: the localized and extended wavefunctions typically do 
overlap giving rise to the repulsion of nearby levels. The formula @ is expected to be valid only 
in the universal regime L <C L*, where L* is a saturation length scale L* ~ -p (d / dk)N (k) 
and I* is the length of the shortest periodic orbit of the desymmetrized billiard, while for L ~ L* 
the shortest periodic orbits of the system start to govern the behavior of a number variance in 

a system-dependent non-universal way ||, ||. Note that L* oc k^ D -^ oc N {D ' 1)/D , where D 
is the number of freedoms, so the universality region of L is much larger in 3D systems than in 
2D systems for comparable sequential quantum numbers: so the quantitative effect of quantum 
eigenstates which are enhanced in the neighborhood of classically invariant sets (scars) on the 
long-range energy level statistics is cleaner in higher dimensions. 

In a conclusion we should stress that using an efficient numerical technique || (see also m) 
we have been able to calculate high-quality energy spectra of a generic 3D chaotic billiard going 
up to 62, 000th excited state within fixed symmetry class. Using sensitive tests on various energy 
level statistics, such as short-range nearest neighbor level spacing distribution and long-range 
number variance, we have detected small but significant deviations from the GOE statistics which 
is an indication of the localization of eigenstates, although the classical dynamics is non-diffusive 
and strongly chaotic. It has been shown that the spectrum can be decomposed with respect to 
long-range correlations to a GOE subsequence of extended chaotic states and an uncorrelated 
(Poisson) subsequence of strongly localized states. The structure of individual eigenfunctions 
a chaotic 3D billiard (|]) in various representations, and examples of strongly localized states, 
are studied in The present study also motivates further study of energy level statistics in 
a generic classically mixed regime, which is also realized by our billiard with smaller absolute 
values of parameters a and b, to test the Berry-Robnik surmise in 3D, and to search for quantal 
consequences of genuine 3D aspects of dynamics, such as Arnold's diffusion and its effect on 
quantum energy level statistics. 
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Figure captions 



Figure 1: The shape of the boundary of the chaotic 3D billiard (|l]) for a = —1/5, b = —12/5. 

Figure 2: The cumulative distributions of (finite-time) maximal liapunov exponents are plotted 
for 7 different lengths of orbits (with unit velocity) l = t = 400, 400^2, 800, 800\/2, 1600, 1600\/2, 3200. 
Each distribution was generated using 10000 orbits with initial conditions chosen randomly and 
uniformly over the phase-space. All distributions have almost the same average maximal lia- 
punov exponent (Xmax) ~ 0.540 with small and decreasing dispersions {(X ma x — (A ma z)) 2 ) = 
0.0567,0.0493,0.0437,0.0391,0.0352,0.0319,0.0290, respectively. The region of small liapunov 
exponents, indicated with a dotted line, is magnified in the inset. The size of a jump of a cu- 
mulative distribution of long-time liapunov exponents for small values of the exponent (which is 
around one tick of the inset = 10~ 3 ) is an estimate (or an upper bound) of the relative volume 
of the regular part of phase space pi . 

Figure 3: The oscillatory part of the level counting function, M(k) —M S mooth(k), is shown for 
the sample I. Local average of the noisy data over 1000 neighboring levels is shown as the small 
band whose width is an estimated statistical error. The fact that the average band overlaps 
with or slightly fluctuates (oscillations are due to shortest periodic orbits) around abscissa on a 
scale of few percents of a level is a very good 'experimental' test of 3D Weyl formula (^). 

Figure 4: The quantum length spectrum D(x) (||) is shown for sample I (a) and sample II (b). 
Vertical dotted lines indicate the lengths of least unstable (with y/\ det M — 1| < 100, where M 
is a monodromy matrix) and shortest (up to 12 bounces) periodic orbits of a desymmetrized 
billiard which match with peaks of quantum length spectrum of the sample I. For the sample I 
(a) we have better resolution that for the sample II (b). The pattern of dots denotes the type of 

stability: . . . elliptic orbit (two pairs of unimodular eigenvalues of M), hyperbolic orbit 

(one pair of unimodular and one pair of real eigenvalues of M), doubly hyperbolic orbit 

(two pairs of real eigenvalues of M), and loxodromic orbit (four complex eigenvalues of 

M). 

Figure 5: The cumulative NNLSD W{S) is shown in (a) for sample I (thin curve) and for 
sample II (thick curve). In the inset we show magnified small S region. The dotted curves are 
the GOE and Poissonian NNLSD. Below (b) we show the deviations from GOE statistics with 
a uniform error and with uniform density of points along abscissa, U(W(S)) — U(Wgoe{S)) 
against W(S), where U{W) = (2/n) arccos \/l ~ W. The thin and thick noisy curve are the 
data for samples I and II, respectively, the dotted and dashed curves are the best fitting Brody 
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and Berry-Robnik distributions, respectively, the thin ones for the sample I and the thick ones 
for the sample II. The horizontal dotted lines indicate ± one-sigma band of estimated statistical 
error for the numerical data. 

Figure 6: We show the number variance S 2 (L) for sample I (thin curve) and sample II (thick 
curve). The dashed curves are the corresponding best fits to Seligman-Verbaarschot formula (||) 
(see text), which has been fitted on the interval 2 < L < 12 for the sample / and on the interval 
2 < L < 40 for the sample II. One should observe excellent agreement with the model (dashed) 
curves (|j) for the spectral sample I up to L ~ 12, and for the spectral sample II up to L ~ 40 
where notable agreement is continued up to L ~ 100. The dotted curve is the number variance 
for GOE. 
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